rm(list=ls())

# load packages 
library(foreign)
library(lmtest)
library(sandwich)
library(stargazer)
library(msm)
library("readxl")
library("ggplot2")
library(gridExtra)
library(readstata13)
library(haven)
library(dplyr)
library(forcats)
library(ggeasy)
library(stargazer)
library(gt)
library(haven)
library(tidyverse)
library(broom)
library(ggplot2)
library(RDHonest)

vindoc_cy_full <- read_csv("Data/vindoc_cy.csv")
chile <- subset(vindoc_cy_full, country_name == "Chile" & year !=1973 & (year >= 1950 & year <= 2000))


## Figure 1

dark_palette <- c("#990000", "#003399", "#006600", "#660066")

inco_plot <- ggplot(chile, aes(x = year)) +
  geom_line(aes(y = v2xed_ed_inpt, color = "Indoctrination Potential"), size = 1, linetype = "solid") +
  geom_line(aes(y = v2xed_ed_cent, color = "Centralization Ed. System"), size = 1, linetype = "solid") +
  geom_line(aes(y = v2xed_ed_inco, color = "Indoctrination Coherence"), size = 1, linetype = "solid") +
  geom_line(aes(y = v2xed_ed_ptcon, color = "Patriotic Indoctrination"), size = 1, linetype = "solid") +
  labs(x = "Year", y = "Indoctrination Scale [0,1]") +
  scale_color_manual(values = dark_palette,
                     labels = c("Indoctrination Potential", "Centralization Ed. System", "Indoctrination Coherence", "Patriotic Indoctrination"),
                     name = NULL) +  # Set name to NULL to remove legend title
  theme_minimal() +
  geom_vline(xintercept = 1972, linetype = "dashed", color = "black") +
  annotate("text", x = 1973, y = Inf, label = "1973", vjust = -0.5, hjust = 0.5, color = "black") +  # Add annotation for the year
  theme(axis.text.x = element_text(size = 14),  # Increase size of x-axis text
        axis.text.y = element_text(size = 14),  # Increase size of y-axis text
        axis.title.x = element_text(size = 14),  # Increase size of x-axis label
        axis.title.y = element_text(size = 14),  # Increase size of y-axis label
        legend.text = element_text(size = 14),   # Increase size of legend text
        legend.title = element_text(size = 14),  # Increase size of legend title
        panel.grid.major = element_line(color = "gray", linetype = "dotted"),
        panel.grid.minor = element_blank()) +  # Add gridlines
  scale_x_continuous(breaks = c(seq(min(chile$year), max(chile$year), by = 15), 1973))  # Add 1973 as a tick mark

ggsave("Output/figure1.png", width = 987, height = 607, units="px", dpi = 72, bg='#ffffff')

print("Figure 1 complete")

## Table 1: Descriptive Statistics

## See 01_results_descriptive.R

## Table 2: RD Estimates 1973 Coup on Ideological Identification

cep_all <- read_dta("Data/cep_all.dta")

# Define dosage samples
dosage_ranges <- c(1, 2, 3, 4, 5)
sample_list <- lapply(dosage_ranges, function(i) subset(cep_all, treatment_dosage_coup >= -i & treatment_dosage_coup <= i))
names(sample_list) <- paste0("sample", dosage_ranges)

# Clustered SE function
cluster_se <- function(model, data, type = "HC1") {
  idx <- as.integer(rownames(model.frame(model)))
  vc  <- sandwich::vcovCL(model, cluster = data$cluster[idx], type = type)
  lmtest::coeftest(model, vcov. = vc)
}

# Generic model runner
run_models <- function(samples, formula, dv_name, table_title, dep_caption, cov_label) {
  models_lm <- lapply(samples, function(s) lm(formula, data = s))
  models_cl <- Map(cluster_se, models_lm, samples)
  sample_sizes <- sapply(models_lm, nobs)
  r_squared <- sapply(models_lm, function(m) summary(m)$r.squared)
  mean_dv <- sapply(models_lm, function(m) mean(m$model[[dv_name]], na.rm = TRUE))
  
  stargazer(models_cl,
            type = "latex",
            title = table_title,
            keep = "^treat1$",
            column.labels = c("54 and 56", "53-57", "52-58", "51-59", "50-60"),
            covariate.labels = cov_label,
            dep.var.caption = dep_caption,
            dep.var.labels = dv_name,
            dep.var.labels.include = FALSE,
            add.lines = list(
              c("Obs.", sample_sizes),
              c("R^2", round(r_squared, 3)),
              c("Mean DV", round(mean_dv, 3))
            ),
            align = TRUE)
}

# ------------------
# Table 2: Left Identification
tab2_left <- run_models(sample_list,
           left ~ treat1 + treatment_dosage_coup + treatment_dosage_coup:treat1 + as.factor(encuesta),
           dv_name = "left",
           table_title = "Table 2: RD Estimates 1973 Coup on Left Identification",
           dep_caption = "Dependent Variable: Left Identification",
           cov_label = "Coup")

# Table 2: Right Identification
tab2_right <- run_models(sample_list,
           right ~ treat1 + treatment_dosage_coup + treatment_dosage_coup:treat1 + as.factor(encuesta),
           dv_name = "right",
           table_title = "Table 2: RD Estimates 1973 Coup on Right Identification",
           dep_caption = "Dependent Variable: Right Identification",
           cov_label = "Coup")

# Table 2: Center Identification
tab2_center <- run_models(sample_list,
           center ~ treat1 + treatment_dosage_coup + treatment_dosage_coup:treat1 + as.factor(encuesta),
           dv_name = "center",
           table_title = "Table 2: RD Estimates 1973 Coup on Center Identification",
           dep_caption = "Dependent Variable: Center Identification",
           cov_label = "Coup")

# Table 2: No Ideology
tab2_noideology <- run_models(sample_list,
           none_pol ~ treat1 + treatment_dosage_coup + treatment_dosage_coup:treat1 + as.factor(encuesta),
           dv_name = "none_pol",
           table_title = "Table 2: RD Estimates 1973 Coup on No Ideological Identification",
           dep_caption = "Dependent Variable: No Ideology",
           cov_label = "Coup")

print("Table 2 complete")

# ------------------
# Table 3: Hard Left
tab3_hardleft <- run_models(sample_list,
           left_hard ~ treat1 + treatment_dosage_coup + treatment_dosage_coup:treat1 + as.factor(encuesta),
           dv_name = "left_hard",
           table_title = "Table 3: RD Estimates 1973 Coup on Hard Left Identification",
           dep_caption = "Dependent Variable: Hard Left Identification",
           cov_label = "Coup")

# Table 3: Soft Left
tab3_softleft <- run_models(sample_list,
           left_soft ~ treat1 + treatment_dosage_coup + treatment_dosage_coup:treat1 + as.factor(encuesta),
           dv_name = "left_soft",
           table_title = "Table 3: RD Estimates 1973 Coup on Soft Left Identification",
           dep_caption = "Dependent Variable: Soft Left Identification",
           cov_label = "Coup")

print("Table 3 complete")

# ------------------
# Table 4: Left Identification by Education
samples_hs <- lapply(sample_list, function(s) subset(s, ed_levels == 4))
samples_college <- lapply(sample_list, function(s) subset(s, ed_levels %in% c(5, 6)))

# High School
tab4_hs <- run_models(samples_hs,
           left ~ treat1 + treatment_dosage_coup + treatment_dosage_coup:treat1 + as.factor(encuesta),
           dv_name = "left",
           table_title = "Table 4: RD Estimates 1973 Coup on Left Identification (High School)",
           dep_caption = "Dependent Variable: Left Identification",
           cov_label = "Coup")

# More than High School
tab4_morehs <- run_models(samples_college,
           left ~ treat1 + treatment_dosage_coup + treatment_dosage_coup:treat1 + as.factor(encuesta),
           dv_name = "left",
           table_title = "Table 4: RD Estimates 1973 Coup on Left Identification (College or More)",
           dep_caption = "Dependent Variable: Left Identification",
           cov_label = "Coup")

print("Table 4 complete")

# ------------------
# Figure 2: Left Identification by Dosage
sample4 <- sample_list[[4]]
sample4$Group <- ifelse(sample4$treatment_dosage_coup <= -1, "Control", ifelse(sample4$treatment_dosage_coup >= 1, "Treatment", "otherwise"))

plot_left <- ggplot(sample4, aes(x = treatment_dosage_coup, y = left, by = Group)) +
  stat_summary(fun = "mean", geom = "point", position = position_dodge(width = 1)) +
  geom_smooth(data = subset(sample4, treatment_dosage_coup <= -1), method = "lm", se = FALSE, color = "grey", linetype = "dashed") +
  geom_smooth(data = subset(sample4, treatment_dosage_coup >= 1), method = "lm", se = FALSE, color = "grey", linetype = "dashed") +
  geom_vline(xintercept = 0, linetype = "dashed", color = "black") +
  labs(x = "Years Exposed to High School Under Pinochet", y = "Left Identification (Average)") +
  theme_minimal() +
  scale_x_continuous(breaks = seq(-4, 4, 1)) +
  theme(axis.title.x = element_text(size = 14),
        axis.title.y = element_text(size = 14),
        axis.text.x = element_text(size = 14),
        axis.text.y = element_text(size = 14))

print(plot_left)
ggsave("Output/figure2.pdf", plot = plot_left, width = 7, height = 5, units = "in")

print("Figure 4 complete")

